#pragma once

#define NMOMMAX 500

float theta[4];

int Nmom;
int Nindexmom;
int *Nmom_per_set;

int ***iP;
double *P2; double *SinP2; double **P; double **SinP; double *SinP4;

int Nmom_eq;
double *P2_eq; double *SinP2_eq; double **P_eq; double **SinP_eq; double *SinP4_eq;

void Momentum(char *momentum_order)
{


  const double PI=3.14159265358979323846;

  int lP[4],size[4],imom=0;
  
  if (strcmp (momentum_order,"TXYZ") == 0)  
  {
	printf("TXYZ momentum order \n");
  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
	 int imom_set=0;

           for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
            for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
             for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
                    for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)

                  {
                    P2[imom]=SinP2[imom]=SinP4[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=T;
                        else size[idir]=L;
                        P[imom][idir]=2*PI*(lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2[imom]+=pow(P[imom][idir],2);
                        SinP[imom][idir]=sin(P[imom][idir]);
                        SinP2[imom]+=pow(SinP[imom][idir],2);
                        SinP4[imom]+=pow(SinP[imom][idir],4);
                      }
                    imom++;
		    imom_set++;
                  }
	   if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
  }
  if (strcmp (momentum_order,"TYXZ") == 0)
  {
        printf("TYXZ momentum order \n");
  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

           for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
            for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
             for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
                    for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)

                  {
                    P2[imom]=SinP2[imom]=SinP4[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=T;
                        else size[idir]=L;
                        P[imom][idir]=2*PI*(lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2[imom]+=pow(P[imom][idir],2);
                        SinP[imom][idir]=sin(P[imom][idir]);
                        SinP2[imom]+=pow(SinP[imom][idir],2);
                        SinP4[imom]+=pow(SinP[imom][idir],4);
                      }
                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
  }

  if (strcmp (momentum_order,"TYZX") == 0)
  {
        printf("TYZX momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
       for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
                  {
                    P2[imom]=SinP2[imom]=SinP4[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=T;
                        else size[idir]=L;
                        P[imom][idir]=2*PI*(lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2[imom]+=pow(P[imom][idir],2);
                        SinP[imom][idir]=sin(P[imom][idir]);
                        SinP2[imom]+=pow(SinP[imom][idir],2);
                        SinP4[imom]+=pow(SinP[imom][idir],4);
                      }

                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
      }
  }
}
void Momentum_Magro(char *momentum_order, int set_p[4])
{

  int lP[4],size[4],imom=0;

  if (strcmp (momentum_order,"TXYZ") == 0)
  {
        printf("TXYZ momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

           for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
            for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
             for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
                    for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)

                  {
                    P2[imom]=SinP2[imom]=SinP4[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=T;
                        else size[idir]=L;
                        P[imom][idir]=2*PI*(set_p[idir]*lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2[imom]+=pow(P[imom][idir],2);
                        SinP[imom][idir]=sin(P[imom][idir]);
                        SinP2[imom]+=pow(SinP[imom][idir],2);
                        SinP4[imom]+=pow(SinP[imom][idir],4);

                      }
                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
 }
  if (strcmp (momentum_order,"TYXZ") == 0)
  {
        printf("TYXZ momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
       for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
        {
                    P2[imom]=SinP2[imom]=SinP4[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=T;
                        else size[idir]=L;
                        P[imom][idir]=2*PI*(set_p[idir]*lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2[imom]+=pow(P[imom][idir],2);
                        SinP[imom][idir]=sin(P[imom][idir]);
                        SinP2[imom]+=pow(SinP[imom][idir],2);
                        SinP4[imom]+=pow(SinP[imom][idir],4);


                      }

                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
  }

  if (strcmp (momentum_order,"TYZX") == 0)
  {
        printf("TYZX momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
       for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
        {
                    P2[imom]=SinP2[imom]=SinP4[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=T;
                        else size[idir]=L;
                        P[imom][idir]=2*PI*(set_p[idir]*lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2[imom]+=pow(P[imom][idir],2);
                        SinP[imom][idir]=sin(P[imom][idir]);
                        SinP2[imom]+=pow(SinP[imom][idir],2);
                        SinP4[imom]+=pow(SinP[imom][idir],4);


                      }

                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
 }
}



void compute_nn(char *momentum_order, double *nn)
{
  int lP[4];
  int imom=0;
  if (strcmp (momentum_order,"TXYZ") == 0)
  {
 for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
     for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
       for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
         for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)

        {
                nn[imom]=4.0;
                if (lP[1]==0) nn[imom]=nn[imom]-1.0;
                if (lP[2]==0) nn[imom]=nn[imom]-1.0;
                if (lP[3]==0) nn[imom]=nn[imom]-1.0;
                imom++;
        }  
   }
  if (strcmp (momentum_order,"TYXZ") == 0)
  {

 for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
       for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
        {
                nn[imom]=4.0;
                if (lP[1]==0) nn[imom]=nn[imom]-1.0;
                if (lP[2]==0) nn[imom]=nn[imom]-1.0;
                if (lP[3]==0) nn[imom]=nn[imom]-1.0;
                imom++;
  }     }

  if (strcmp (momentum_order,"TYZX") == 0)
  {

 for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
       for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
	{
		nn[imom]=4.0;
		if (lP[1]==0) nn[imom]=nn[imom]-1.0;
                if (lP[2]==0) nn[imom]=nn[imom]-1.0;
                if (lP[3]==0) nn[imom]=nn[imom]-1.0;
		imom++;
  }	}
}


int  Number_of_Equivalent_Momenta(double *Eq_for_P2)
{
int eq_imom_find;
  int imom_eq=0;
  for (int imom=0; imom<Nmom; imom++)
  {
	eq_imom_find=0;
        for (int imom_loop=0; imom_loop<imom; imom_loop++)
        {
                if (fabs(Eq_for_P2[imom_loop]-Eq_for_P2[imom])<0.000001)
                {
                         //printf("imom = %d has the same SinP2 as imom=%d\n",imom,imom_loop);
                         eq_imom_find=1;
                 }
         }
         // Si no hemos encontrado ningun momento equivalente, este debe ser uno nuevo
         if (eq_imom_find==0)   imom_eq++; 
   }
   int  Nmom_eq=imom_eq;

          P2_eq=(double*)malloc(sizeof(double)*(Nmom_eq));
          SinP2_eq=(double*)malloc(sizeof(double)*(Nmom_eq));
          P_eq=(double**)malloc(sizeof(double*)*(Nmom_eq));
          for (int ip=0 ; ip<Nmom_eq; ip++) P_eq[ip]=(double*)malloc(sizeof(double)*(4));
          SinP_eq=(double**)malloc(sizeof(double*)*(Nmom_eq));
          for (int ip=0 ; ip<Nmom_eq; ip++) SinP_eq[ip]=(double*)malloc(sizeof(double)*(4));
          SinP4_eq=(double*)malloc(sizeof(double)*(Nmom_eq));
   return Nmom_eq;
}

/*void  Equivalent_Momenta(double *Eq_for_P2)
{
int first_eq_imom;
int eq_imom_find=0;
  int imom_eq=0;
  for (int imom=0; imom<Nmom; imom++)
  {
	eq_imom_find=0;
        for (int imom_loop=0; imom_loop<imom; imom_loop++)
        {
                if (fabs(Eq_for_P2[imom_loop]-Eq_for_P2[imom])<0.000001)
                {
			  //printf("imom=%d has the same p2 as imom=%d\n",imom, imom_loop);
                          if (eq_imom_find==0) first_eq_imom=imom_loop; 
                          eq_imom_find=1;
                 }
         }
	 // Si este momento no se habia contabilizado aun guardarlo en imom_eq
            if (eq_imom_find==0)
	    {
		// printf("imom=%d  with p2=%0.5f is new saved as imom=%d\n",imom,Eq_for_P2[imom],imom_eq);
		 P2_eq[imom_eq]=P2[imom];
		 SinP2_eq[imom_eq]=SinP2[imom];
		 for (int idir=0; idir<4; idir++)  P_eq[imom_eq][idir]=P[imom][idir];
		 for (int idir=0; idir<4; idir++)  SinP_eq[imom_eq][idir]=SinP[imom][idir];
		 SinP4_eq[imom_eq]=SinP4[imom];
	    }
         // Si no hemos encontrado ningun momento equivalente, este debe ser uno nuevo
         if (eq_imom_find==0) imom_eq++; 
	 // Si hemos encontrado alguno equivalente (eq_imom_find==1) entonces este no cuenta
   }

}
*/

void  Average_over_Equivalent_Momenta(double *Eq_for_P2, double *Z, double *Z_eq)
{
int first_eq_imom;
int eq_imom_find=0;

// Defino estas variables auxiliares con un numero maximo de momentos
// Con punteros me da error

int number_eq[NMOMMAX];
int list_eq_imom[NMOMMAX];
double ZAUX_eq[NMOMMAX];
//

int imom_eq=0;
  for (int imom=0; imom<Nmom; imom++)
  {
	eq_imom_find=0;
	for (int imom_loop=0; imom_loop<imom; imom_loop++)
        {
                if (fabs(Eq_for_P2[imom_loop]-Eq_for_P2[imom])<0.000001)
                {
                          if (eq_imom_find==0) first_eq_imom=imom_loop; // Debe sumarse al primero
                          eq_imom_find=1;
			  //printf("imom=%d is equivalent to %d\n",imom,first_eq_imom);
                 }
         }
	 // Si no hemos encontrado ningun momento equivalente el valor de Z[imom] se guarda en un nuevo momento
         // Si no hemos encontrado ningun momento equivalente, este debe ser uno nuevo por lo tanto inicializo el number_eq para ese momento
        // Si lo hemos encontrado sumo 1 al number_eq para ese imom_eq

         if (eq_imom_find==0)
	 {
		 ZAUX_eq[imom_eq]=Z[imom];
		 number_eq[imom_eq]=1; 
		 list_eq_imom[imom]=imom_eq;
//		 printf("imom %d is new,Z=%8.5f=%8.5f saving in imom_eq=%d=%d\n",imom,Z[imom],ZAUX_eq[imom_eq],imom_eq,list_eq_imom[imom]);
		 imom_eq++;
          }
         if (eq_imom_find==1)
	 {
		ZAUX_eq[list_eq_imom[first_eq_imom]]=ZAUX_eq[list_eq_imom[first_eq_imom]]+Z[imom];
		number_eq[list_eq_imom[first_eq_imom]]=number_eq[list_eq_imom[first_eq_imom]]+1;
//		printf("imom=%d is equivalent to %d corresponding to imom_eq=%d, Z=%0.5f  with number_eq %d\n",imom,first_eq_imom, list_eq_imom[first_eq_imom],ZAUX_eq[list_eq_imom[first_eq_imom]],number_eq[list_eq_imom[first_eq_imom]]);
	 }


   }
   for (int imom_eq=0; imom_eq<Nmom_eq; imom_eq++){
         Z_eq[imom_eq]=ZAUX_eq[imom_eq]/number_eq[imom_eq];

    }
}
void  Average_over_Equivalent_Momenta_LT(double *Eq_for_P2, double *Z, double *Z_eq, int Nmom_eq_X)
{
int first_eq_imom;
int eq_imom_find=0;

// Defino estas variables auxiliares con un numero maximo de momentos
// Con punteros me da error

int number_eq[NMOMMAX];
int list_eq_imom[NMOMMAX];
double ZAUX_eq[NMOMMAX];
//

int imom_eq=0;
  for (int imom=0; imom<Nmom; imom++)
  {
        eq_imom_find=0;
        for (int imom_loop=0; imom_loop<imom; imom_loop++)
        {
                if (fabs(Eq_for_P2[imom_loop]-Eq_for_P2[imom])<0.000001)
                {
                          if (eq_imom_find==0) first_eq_imom=imom_loop; // Debe sumarse al primero
                          eq_imom_find=1;
                          //printf("imom=%d is equivalent to %d\n",imom,first_eq_imom);
                 }
         }
         // Si no hemos encontrado ningun momento equivalente el valor de Z[imom] se guarda en un nuevo momento
         // Si no hemos encontrado ningun momento equivalente, este debe ser uno nuevo por lo tanto inicializo el number_eq para ese momento
        // Si lo hemos encontrado sumo 1 al number_eq para ese imom_eq

         if (eq_imom_find==0)
         {
                 ZAUX_eq[imom_eq]=Z[imom];
                 number_eq[imom_eq]=1;
                 list_eq_imom[imom]=imom_eq;
//               printf("imom %d is new,Z=%8.5f=%8.5f saving in imom_eq=%d=%d\n",imom,Z[imom],ZAUX_eq[imom_eq],imom_eq,list_eq_imom[imom]);
                 imom_eq++;
          }
         if (eq_imom_find==1)
         {
                ZAUX_eq[list_eq_imom[first_eq_imom]]=ZAUX_eq[list_eq_imom[first_eq_imom]]+Z[imom];
                number_eq[list_eq_imom[first_eq_imom]]=number_eq[list_eq_imom[first_eq_imom]]+1;
//              printf("imom=%d is equivalent to %d corresponding to imom_eq=%d, Z=%0.5f  with number_eq %d\n",imom,first_eq_imom, list_eq_imom[first_eq_imom],ZAUX_eq[list_eq_imom[first_eq_imom]],number_eq[list_eq_imom[first_eq_imom]]);
         }


   }
   for (int imom_eq=0; imom_eq<Nmom_eq_X; imom_eq++){
         Z_eq[imom_eq]=ZAUX_eq[imom_eq]/number_eq[imom_eq];

    }
}

void  Average_over_Equivalent_Momenta_c(double *Eq_for_P2, complex *Z, complex *Z_eq)
{
int first_eq_imom;
int eq_imom_find=0;

// Defino estas variables auxiliares con un numero maximo de momentos
// Con punteros me da error

int number_eq[NMOMMAX];
int list_eq_imom[NMOMMAX];
complex ZAUX_eq[NMOMMAX];
//

int imom_eq=0;
  for (int imom=0; imom<Nmom; imom++)
  {
        eq_imom_find=0;
        for (int imom_loop=0; imom_loop<imom; imom_loop++)
        {
                if (fabs(Eq_for_P2[imom_loop]-Eq_for_P2[imom])<0.000001)
                {
                          if (eq_imom_find==0) first_eq_imom=imom_loop; // Debe sumarse al primero
                          eq_imom_find=1;
                          //printf("imom=%d is equivalent to %d\n",imom,first_eq_imom);
                 }
         }
         // Si no hemos encontrado ningun momento equivalente el valor de Z[imom] se guarda en un nuevo momento
         // Si no hemos encontrado ningun momento equivalente, este debe ser uno nuevo por lo tanto inicializo el number_eq para ese momento
        // Si lo hemos encontrado sumo 1 al number_eq para ese imom_eq

         if (eq_imom_find==0)
         {
                 ZAUX_eq[imom_eq]=Z[imom];
                 number_eq[imom_eq]=1;
                 list_eq_imom[imom]=imom_eq;
                 //printf("imom %d is new,Z=%8.5f=%8.5f saving in imom_eq=%d=%d\n",imom,Z[imom],ZAUX_eq[imom_eq],imom_eq,list_eq_imom[imom]);
                 imom_eq++;
          }
         if (eq_imom_find==1)
         {
                ZAUX_eq[list_eq_imom[first_eq_imom]]=ZAUX_eq[list_eq_imom[first_eq_imom]]+Z[imom];
                number_eq[list_eq_imom[first_eq_imom]]=number_eq[list_eq_imom[first_eq_imom]]+1;
                //printf("imom=%d is equivalent to %d corresponding to imom_eq=%d, Z=%0.5f  with number_eq %d\n",imom,first_eq_imom, list_eq_imom[first_eq_imom],ZAUX_eq[list_eq_imom[first_eq_imom]],number_eq[list_eq_imom[first_eq_imom]]);
         }


   }
   for (int imom_eq=0; imom_eq<Nmom_eq; imom_eq++){
         Z_eq[imom_eq]=ZAUX_eq[imom_eq]/number_eq[imom_eq];

    }
}

void read_momenta_data(char *file)
{

           char input_momdat[1024];
           sprintf(input_momdat,file);
		printf("Read %s\n",file);
           FILE *momdat=fopen(input_momdat,"r");

           fscanf(momdat,"%f %f %f %f",
           &theta[0],&theta[1],
           &theta[2],&theta[3]
           );
           fscanf(momdat,"%i",&Nindexmom);

	  iP=(int***)malloc(sizeof(int**)*(Nindexmom));
	  for(int imom=0;imom<Nindexmom;imom++) iP[imom]=(int**)malloc(sizeof(int*)*(4));
  	  for(int imom=0;imom<Nindexmom;imom++) for (int idir=0; idir<4; idir++) iP[imom][idir]=(int*)malloc(sizeof(int)*(2));


           for(int inm=0; inm<Nindexmom; inm++)
            {
              fscanf(momdat,"%i %i %i %i %i %i %i %i",
                   &iP[inm][0][0],&iP[inm][0][1],
                   &iP[inm][1][0],&iP[inm][1][1],
                   &iP[inm][2][0],&iP[inm][2][1],
                   &iP[inm][3][0],&iP[inm][3][1]
                 );
            }

          fclose(momdat);

	  Nmom_per_set=(int*)malloc(sizeof(int)*(Nindexmom));
          for(int imom=0;imom<Nindexmom;imom++){
                        Nmom_per_set[imom]=1;
                        for(int idir=0; idir<4; idir++) Nmom_per_set[imom]=Nmom_per_set[imom]*(iP[imom][idir][1]-iP[imom][idir][0]+1);
             }
           int Nmom_tot=0;
           for (int imom=0; imom<Nindexmom;imom++) Nmom_tot=Nmom_tot+Nmom_per_set[imom];

           Nmom=Nmom_tot;
          // Momentum variables definition
          P2=(double*)malloc(sizeof(double)*(Nmom));
          SinP2=(double*)malloc(sizeof(double)*(Nmom));
          P=(double**)malloc(sizeof(double*)*(Nmom));
          for (int imom=0 ; imom<Nmom; imom++) P[imom]=(double*)malloc(sizeof(double)*(4));
          SinP=(double**)malloc(sizeof(double*)*(Nmom));
          for (int imom=0 ; imom<Nmom; imom++) SinP[imom]=(double*)malloc(sizeof(double)*(4));
          SinP4=(double*)malloc(sizeof(double)*(Nmom));
          

}
void Momentum_LT(double *P2_X, double *SinP2_X, double **P_X, double **SinP_X, double *SinP4_X, char *momentum_order, int LX, int TX)
{


  const double PI=3.14159265358979323846;

  int lP[4],size[4],imom=0;
  
  if (strcmp (momentum_order,"TXYZ") == 0)  
  {
	printf("TXYZ momentum order \n");
  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
	 int imom_set=0;

           for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
            for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
             for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
                    for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)

                  {
                    P2_X[imom]=SinP2_X[imom]=SinP4_X[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=TX;
                        else size[idir]=LX;
                        P_X[imom][idir]=2*PI*(lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2_X[imom]+=pow(P_X[imom][idir],2);
                        SinP_X[imom][idir]=sin(P_X[imom][idir]);
                        SinP2_X[imom]+=pow(SinP_X[imom][idir],2);
                        SinP4_X[imom]+=pow(SinP_X[imom][idir],4);
                      }
                    imom++;
		    imom_set++;
                  }
	   if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
  }
  if (strcmp (momentum_order,"TYXZ") == 0)
  {
        printf("TYXZ momentum order \n");
  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

           for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
            for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
             for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
                    for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)

                  {
                    P2_X[imom]=SinP2_X[imom]=SinP4_X[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=TX;
                        else size[idir]=LX;
                        P_X[imom][idir]=2*PI*(lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2_X[imom]+=pow(P_X[imom][idir],2);
                        SinP_X[imom][idir]=sin(P_X[imom][idir]);
                        SinP2_X[imom]+=pow(SinP_X[imom][idir],2);
                        SinP4_X[imom]+=pow(SinP_X[imom][idir],4);
                      }
                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
  }

  if (strcmp (momentum_order,"TYZX") == 0)
  {
        printf("TYZX momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
       for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
                  {
                    P2_X[imom]=SinP2_X[imom]=SinP4_X[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=TX;
                        else size[idir]=LX;
                        P_X[imom][idir]=2*PI*(lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2_X[imom]+=pow(P_X[imom][idir],2);
                        SinP_X[imom][idir]=sin(P_X[imom][idir]);
                        SinP2_X[imom]+=pow(SinP_X[imom][idir],2);
                        SinP4_X[imom]+=pow(SinP_X[imom][idir],4);
                      }

                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
      }
  }
}

void Momentum_Magro_LT(double *P2_X, double *SinP2_X, double **P_X, double **SinP_X, double *SinP4_X, char *momentum_order, int set_p[4], int LX, int TX)
{

  int lP[4],size[4],imom=0;

  if (strcmp (momentum_order,"TXYZ") == 0)
  {
        printf("TXYZ momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

           for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
            for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
             for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
                    for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)

                  {
                    P2_X[imom]=SinP2_X[imom]=SinP4_X[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=TX;
                        else size[idir]=LX;
                        P_X[imom][idir]=2*PI*(set_p[idir]*lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2_X[imom]+=pow(P_X[imom][idir],2);
                        SinP_X[imom][idir]=sin(P_X[imom][idir]);
                        SinP2_X[imom]+=pow(SinP_X[imom][idir],2);
                        SinP4_X[imom]+=pow(SinP_X[imom][idir],4);

                      }
                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
 }
  if (strcmp (momentum_order,"TYXZ") == 0)
  {
        printf("TYXZ momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;

  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
       for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
        {
                    P2_X[imom]=SinP2_X[imom]=SinP4_X[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=TX;
                        else size[idir]=LX;
                        P_X[imom][idir]=2*PI*(set_p[idir]*lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2_X[imom]+=pow(P_X[imom][idir],2);
                        SinP_X[imom][idir]=sin(P_X[imom][idir]);
                        SinP2_X[imom]+=pow(SinP_X[imom][idir],2);
                        SinP4_X[imom]+=pow(SinP_X[imom][idir],4);


                      }

                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
  }

  if (strcmp (momentum_order,"TYZX") == 0)
  {
        printf("TYZX momentum order \n");

  for (int indexmom=0; indexmom<Nindexmom; indexmom++)
  {
         int imom_set=0;
  for(lP[0]=iP[indexmom][0][0];lP[0]<=iP[indexmom][0][1];lP[0]++)
    for(lP[2]=iP[indexmom][2][0];lP[2]<=iP[indexmom][2][1];lP[2]++)
     for(lP[3]=iP[indexmom][3][0];lP[3]<=iP[indexmom][3][1];lP[3]++)
       for(lP[1]=iP[indexmom][1][0];lP[1]<=iP[indexmom][1][1];lP[1]++)
        {
                    P2_X[imom]=SinP2_X[imom]=SinP4_X[imom]=0;
                    for(int idir=0;idir<4;idir++)
                      {
                        if (idir==0) size[idir]=TX;
                        else size[idir]=LX;
                        P_X[imom][idir]=2*PI*(set_p[idir]*lP[idir]+theta[idir]*0.5)/size[idir];
                        //P[idir][imom]=2*PI*lP[idir]/size[idir] - PI/2; // checking Marianne
                        P2_X[imom]+=pow(P_X[imom][idir],2);
                        SinP_X[imom][idir]=sin(P_X[imom][idir]);
                        SinP2_X[imom]+=pow(SinP_X[imom][idir],2);
                        SinP4_X[imom]+=pow(SinP_X[imom][idir],4);


                      }

                    imom++;
                    imom_set++;
                  }
           if (imom_set != Nmom_per_set[indexmom]) printf("Error in momentum routine\n");
  }
 }
}


